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We report a variational calculation of ground state energies and radii of 4 Hejv droplets (3 < N < 40), 
using the atom-atom interaction HFD-B(HE). The trial wave function has a simple structure, 
combining two- and three-body correlation functions coming from a translationally invariant 
configuration-interaction description, and Jastrow-type short-range correlations. The calculated 
ground state energies differ by around 2% from the diffusion Monte Carlo results. 



PACS numbers: 36.40.-c 61.46. +w 

The research on liquid helium clusters has attracted a 
great interest both experimentally and theoretically jl],|| . 
This research allows for the analysis of the evolution of 
various physical properties for increasing size of the sys- 
tem, going from single atoms to the bulk. Helium clusters 
are expected to remain liquid under all conditions of for- 
mation, offering thus the possibility to study finite-size 
effects in the superfluid state ||. Moreover, it has been 
suggested that Bose condensation could be detected by 
means of helium atom-cluster collisions [Q . The experi- 
mental research has faced the difficulties of size detecting 
the clusters. Recently, molecular beam diffraction from 
a transmission grafting Q has proven to be successful to 
detect even the 4 He dimer ||, giving further impetus to 
the study of helium clusters. 

As the atom-atom interaction is well-known and rela- 
tively simple, the solution of the Schrodinger equation 
has been obtained using several microscopic methods, 
mainly based on Monte Carlo techniques j7|J^,p|, |To| , |TT| | . 
Variational Monte Carlo (VMC) calculations, using a 
Jastrow-like ansatz for the many-body wave function, are 
currently used for systems dominated by strong short- 
range interactions. The VMC wave functions are the 
input for diffusion Monte Carlo (DMC), Green function 
Monte Carlo or path integral techniques which provide 
essentially exact, within statistical errors, ground state 
energies of 4 He clusters at zero temperature. Conversely, 
these calculations constitute a useful benchmark to test 
other many-body methods. 

In this work we present a new method of obtaining 
high-quality variational wave functions to describe small 
bosonic clusters. The basic idea is to write the trial wave 
function as the product of three terms, each of them 
with well defined roles. The first term is the familiar 
two-body Jastrow correlation factor, which controls the 
strong atom-atom repulsion at very short distances. The 
second term is related to a single-particle description of 
the cluster and is written as the product of N single 
particle wave functions referred to the center-of-mass co- 



ordinate; it provides the required confinement of the con- 
stituents and fixes basically the size of the droplet. Fi- 
nally, the third term corresponds to a special version of 
the configuration interaction (CI) expansion describing 
two- and three-particle excitations, its role being to in- 
corporate fine details to the wave function for medium 
and long ranges, as well as some collective effects. In 
summary, the trial wave functions we shall consider to 
describe the ground state of 4 He at clusters is 



*(R)=n<?(^)$(R) 



z<3 



i<j i<j<k 

where R represents all single-particle coordinates. The 
two-particle function g(rij) is the two-body Jastrow cor- 
relation. We shall use the form 



,,(r) =o_x.»( -i f£ 



(2) 



introduced by McMillan [Q many years ago for the de- 
scription of the homogeneous liquid using a 12-6 Lennard- 
Jones interaction, with the value v — 5 as required by 
Kato's cusp condition. Many VMC calculations have 
shown that this form is also appropriate when dealing 
with interactions of the Aziz type. For instance, in Ref. 
]l0| fixed values b = 3A, v = 5 have been used for the 
interaction HFDHE2 independently of the number 
of constituents in the cluster. It makes sense to keep the 
same values for all the clusters, because this pair corre- 
lation should be determined only by the two-body inter- 
action potential. We have slightly modified these values, 
employing b — 2.95Aand v = 5.2 in our calculations with 
the interaction HFD-B(HE) @. Note that, contrarily 
to other studies, we are considering for the Jastrow term 
only the dominant part at very short-ranges. 
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The function $(R) is the product of N single-particle 
wave functions 0(r, — Rcm), referred to the center-of- 
mass of the system. If is a gaussian (in other words, 
the Is state of a harmonic oscillator potential) , it may be 
written as 



$(R)=nex P - — 



i<3 



(3) 



depending on a size parameter a. Because of its Jastrow- 
like form we may absorb it into the definition of the Jas- 
trow correlation factor, and define 



(4) 



Our contribution to the variational description of the 
clusters is the term within parentheses in Eq.M), with the 
pair fa and triplet f% correlations, which we have called 
above configuration interaction and will be abridged as 
CI2 and CI3 respectively. In previous works (l7) we have 
extensively applied the CI2 scheme to the description of 
light nuclear systems, with less singular interactions than 
the usual atom-atom potentials. The method is based on 
a linearized version of the coupled cluster method (CCM) 
flEf at the SUB(2) and SUB(3) truncation approximation 
levels, restricted to translationally invariant excitations 
up to three particles-three holes. The linearized version 
of the CCM has been shown to be equivalent to a special 
configuration interaction scheme fl6[| , much more efficient 
than the usual CI methods when dealing with realistic in- 
teractions. With respect to CCM, the CI scheme loses the 
cluster property and the correct scaling with the number 
of particles, being unappropriate for extended systems. 
However it is much simpler, as far as non-linear terms of 
the CCM expansion are absent. 

The present mixed J-CI scheme may be compared with 
the correlated basis functions (CBF) method as devel- 
oped by Feenberg, Clark and Krotscheck |Q for ex- 
tended systems, in which Jastrow correlations are com- 
bined with non-orthogonal perturbation theory. The key 
concept behind both CBF and J-CI methods is to ap- 
ply a divide-and-conquer strategy to the determination 
of the correlated wave function. It is known that Jastrow 
correlations are very well suited to treat the effects of the 
(strongly repulsive) short-range part of the interaction, 
whereas a CI approach can effectively deal with medium- 
and long-range correlations. In our approach, we sub- 
stitute the standard configuration-interaction component 
by our conceptually equivalent but significantly more ef- 
fective translationally invariant CI scheme. Note that 
we are combining additive CI and multiplicative Jastrow 
correlations. 

The CI correlation functions fa and fa could be deter- 
mined by minimizing the ground state energy, resulting 
in a system of coupled integro-differential equations for 
fa and fa. As a practical alternative we have expanded 
these functions in a set of gaussians: 



1 + E h ( r « ) + E h ( r v ' nk ' r 3 k ) = 

i<j i<j<k 

E °p^S J £ e-^e-^e-^* 1 , (5) 

p<q<r— 1 I i<j<k J 

where S indicates symmetrization with respect to the 
particle labels. The gaussian expansion has proven to 
be a very accurate representation of the correlation func- 
tions when both negative and positive values of j3 are 
included, with the only restriction of having a square in- 
tegrable wave function. Among the set of parameters 
{(3p\ we fix one of them, say /3±, to zero. In this way we 
cover the three possibilities contained in the ansatz wave 
function: restricting the three labels p,q,r to 1, we deal 
(up to a normalization constant) with the correlated state 
$j; keeping only two of these labels equal to 1, amounts 
to put on top of that state the linear two-body correla- 
tions; finally the unrestricted choice corresponds to the 
complete ansatz. Calculations presented later on will be 
referred as J, J-CI2 and J-CI3 results, respectively. 

A convenient short-hand notation for the gaussian ex- 
pansion is the following: 



(6) 



where the subindex /i refers to the ordered set (p < q < 
r), and 

F M (R) = S J E ex P " Ar4 " Prr%) \ (7) 

[ i<j<k J 

The variational determination of the energy reduces to 
the solution of a generalized eigenvalue problem, which 
can be stated as follows: 

E (^"Mi ,M2 + H*i .Ma ) = E E ^ ^ ^ 2 

fJ-2 M2 

The matrix elements of the norm and the potential en- 
ergy are the integrals 

N^ 2 = J dRlMR^tRJF^CR) (9) 

and 

V^ 2 = f dR|$,/(R)| 2 F* i (R) E V(r mn )F, 2 (R) 

(10) 

where V(r) is the two-body interaction potential. For the 
matrix elements of the kinetic energy operator we choose 
to write them as 



= J dR|$,/(R)| 2 F* i (R) 
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since we shall use the positive definite function |$j(R)| 2 
as the guide of a Metropolis random walk fl9fl . Note 
that no substraction of the center-of-mass contribution 
is necessary, as we are using a translationally invariant 
wave function. The number of unkown amplitudes C M 
is ( Ar ' 3 3 +2 ). As mentioned above, the J and J-CI2 cases 
can be extracted easily from the general wave function: 
indeed, they correspond to the first d x d block of the 
matrices and the first d amplitudes, with d — 1 and d = 
N/3, respectively. 

The last point to be discussed is the selection of the 
set of parameters {f3 p }. Our previous experience with nu- 
clear systems indicates that this selection is not very crit- 
ical, as far as a sufficiently large interval of length ranges 
{l/yf/3 p ) is included. In the calculations presented be- 
low we have used the set {(3 p /a 2 } = {0, —0.05,0.5, 1,4}. 
Large values of (i p should not be used, to avoid competi- 
tion with the Jastrow factor. 

The computational algorithm consists in carrying out 
a random walk guided by |<I>j(R)| 2 as the probability dis- 
tribution function to evaluate the norm, potential and ki- 
netic matrices as given in eqs.(p[-pT|). The only adjustable 
parameter to be determined by minimization is a, since 
the amplitudes C M are self-adjustable. A final warning 
is in order. It is well known that calculations involv- 
ing the solution of a generalized eigenvalue problem like 
eq.(|^) are prone to numerical instabilities, because the 
over-complete basis giving rise to the overlap matrix Af 
may have a determinant close to zero. In our case, the 
Monte Carlo evaluation of the overlap matrix ([)[) results 
always in a positive matrix, and the possible numerical 
instabilities are related only to numerical rounding er- 
rors. 

In Table I are collected our result s for bin ding energies 
and unit radii, defined as ro = yfbmr^/N 1 ^, where 
(r 2 ) is the square mean radius. Columns labelled J, J- 
CI2 and J-CI3 correspond to the sequence of trial wave 
functions previously explained. Columns labelled VMC 
and DMC display the results obtained by other authors 
using the same interaction. 

Our lowest order approximation (J) should be com- 
pared with column VMC. It can be seen that our crude 
choice for the two-body Jastrow correlation (^J) produces 
less binding than the quoted VMC results. The reason 
of these differences is quite clear: once the parameters 
b and v have been fixed, there remains a single varia- 
tional parameter (a) in <I>j. In contrast, the referred 
VMC calculations use more elaborate forms for the trial 
wave function, including three-body correlation functions 
in some cases, and containing from 5 to 10 parameters. 
Note however, that these differences decrease when N 
increases. 

The inclusion of the CI correlation function /2 changes 
substantially the results. It can be seen that figures in 
columns VMC and J-CI2 are basically the same, within 
statistical errors, except for the cluster with N = 40. 
The J-CI2 values have been obtained by minimizing the 



energy with respect to a and the set of five amplitudes 
describing /2. The resulting values of a are given in Ta- 
ble I. In some cases we have also performed calculations 
including up to nine gaussians. It turns out that the 
binding energy does not change significantly, in spite of 
a considerably increase of the computing time. 

We finally refer to our best variational results, namely 
J-CI3. There are no adjustable parameters in J-CI3, be- 
cause we took the value of a which minimizes the J- 
CI2 energy and the same set of five (3 exponents, the 
amplitudes being determined by solving the generalized 
eigenvalue problem. It can be seen in Table I that the 
ground state energies provided by the J-CI3 scheme are 
above DMC ones by about 2%, except at the extremes 
of the table. For N = 3 both J-CI3 and DMC ener- 
gies agree within statistical errors, and for N — 20 and 
40 the differences are at the level of 4%. In any case, 
the J-CI3 description always shows a large improvement 
with respect to the more elaborate VMC calculations dis- 
played in Table I. Moreover, the values obtained for the 
unit radius are very close to the values given by DMC 
calculations. 

It is interesting to show the sensitivity to the value of 
the size parameter a of the three levels of approximation. 
As a typical example, in Figure 1 are shown our results 
for the cluster with N = 8. For completeness this figure 
includes also two horizontal dotted lines, which represent 
the optimal VMC and DMC energies, and which are not 
related in any form to our parameter a. One observes 
that even if there is a strong dependence of J energies 
with a, this dependence almost disappears when / 2 and 
/3 correlations are included. The flatness of J-CI2 and 
J-CI3 energies is an indication of the completeness of the 
basis used to describe these correlations. 

In conclusion, we have presented in this paper a vari- 
ational wave function for 4 He clusters. A simple Jas- 
trow function of McMillan type is sufficient to screen 
the strong short-range atom-atom repulsion, whereas the 
remaining medium- and long-range correlations seem to 
be adequately described by linear CI correlation func- 
tions. These functions are authomatically determined 
once given the size parameter a, which is the only ad- 
justable parameter. We have shown that these vari- 
ational wave functions produce high-quality results for 
small clusters. The method is likely applicable to other 
bosonic systems. 
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FIG. 1. Ground state energy per particle of cluster 4 Hes, as 
a function of the inverse length parameter a, for the sequence 
of trial wave functions J, J-CI2 and J-CI3. VMC and DMC 
results are taken from Ref. fll|. 
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TABLE I. Ground state energies per particle and unit radii of 4 Hejv droplets, using the interaction HFD-B(HE). VMC and 
DMC energies are taken from [IllJ for N < 10, and |2j otherwise. DMC radii are taken from bj. 





E/N (K) 




ro 


(A) 


N 


J 


J-CI2 


J-CI3 


VMC 


DMC 




J-CI3 


DMC 


3 


+.0136(4) 


-.0373(8) 


-.0430(10) 


-.0385(3) 


-.0436(2) 


.18 


5.4(2) 


5.59 


4 


-.0835(7) 


-.1326(13) 


-.1398(15) 


-.1333(10) 


-.1443(2) 


.22 


4.11(11) 


4.13 


5 


-.1984(8) 


-.2502(14) 


-.2616(13) 


-.2506(3) 


-.2670(3) 


.25 


3.66(7) 


3.65 


6 


-.3187(11) 


-.3678(12) 


-.3868(11) 


-.3676(2) 


-.3950(2) 


.27 


3.42(6) 




7 


-.4360(9) 


-.4823(13) 


-.5081(12) 


-.4845(4) 


-.5206(4) 


.27 


3.31(5) 


3.22 


8 


-.5489(11) 


-.5949(14) 


-.6289(13) 


-.5919(5) 


-.6417(4) 


.28 


3.18(4) 




9 


-.6531(10) 


-.6995(13) 


-.7392(12) 


-.6924(5) 


-.7563(6) 


.28 


3.11(3) 




10 


-.7522(11) 


-.7945(16) 


-.8484(19) 


-.7916(7) 


-.8654(7) 


.30 


3.01(3) 




14 


-1.0908(12) 


-1.130(2) 


-1.215(2) 


-1.1290(7) 


-1.2478(12) 


.28 


2.91(2) 


2.83 


20 


-1.4743(16) 


-1.5106(2) 


-1.6336(15) 


-1.510(2) 


-1.688(2) 


.30 


2.727(14) 


2.69 


40 


-2.218(3) 


-2.273(3) 


-2.4563(14) 


-2.430(2) 


-2.575(3) 


.28 


2.578(8) 
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